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Abstract 



OO 

ly-j , This Mathematica 5.2 package 1 is a simulation of a Quantum Computer. The pro- 

gram provides a modular, instructive approach for generating the basic elements 
that make up a quantum circuit. The main emphasis is on using the density ma- 
S^h. trix, although an approach using state vectors is also implemented in the package. 

The package commands are defined in Q density. m which contains the tools needed 
in quantum circuits, e.g. multiqubit kets, projectors, gates, etc. Selected examples 
of the basic commands are presented here and a tutorial notebook, Tutorial. nb is 
provided with the package (available on our website) that serves as a full guide to 
the package. Finally, application is made to a variety of relevant cases, including 
^ ■ Teleportation, Quantum Fourier transform, Grover's search and Shor's algorithm, in 

separate notebooks: QFT.nb, Teleportation. nb, Grover.nb and Shor.nb where each 
algorithm is explained in detail. Finally, two examples of the construction and ma- 
nipulation of cluster states, which are part of "one way computing" ideas, are in- 
cluded as an additional tool in the notebook Cluster. nb. A Mathematica palette 
containing most commands in QDENSITY is also included: QDENSpalette.nb . 



1 QDENSITY is available at http://www.pitt.edu/~tabakin/QDENSITY 
Preprint submitted to Elsevier Science 1 February 2008 



Program Summary 



Title of program: QDENSITY 
Catalogue identifier: 

Program summary URL: |http: / / cpc.cs.qub.ac.uk/summaries 

Program available from: CPC Program Library, Queen's University of Belfast, 

N. Ireland 

Operating systems: Any which supports Mathematica; tested under Microsoft 
Windows XP, Macintosh OS X, and Linux FC4. 
Programming language used: Mathematica 5.2 

Number of bytes in distributed program, including test code and documenta- 
tion: 1185886 
Distribution format: tar.gz 

Nature of Problem: Analysis and design of quantum circuits, quantum algo- 
rithms and quantum clusters. 

Method of Solution: A Mathematica package is provided which contains com- 
mands to create and analyze quantum circuits. Several Mathematica note- 
books containing relevant examples: Teleportation, Shor's Algorithm and Grover's 
search are explained in detail. A tutorial, Tutorial. nb is also enclosed. 
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1 INTRODUCTION 

There is already a rich Quantum Computing (QC) literature [1] which holds 
forth the promise of using quantum interference and superposition to solve 
otherwise intractable problems. The field has reached the point that experi- 
mental realizations are of paramount importance and theoretical tools towards 
that goal are needed: to gauge the efficacy of various approaches, to under- 
stand the construction and efficiency of the basic logical gates, and to delineate 
and control environmental decoherence effects. 

In this paper, a Mathematica [2] package provides a simulation of a Quantum 
Computer that is both flexible and an improvement over earlier such works [3] . 
It is a bona fide simulation in that its success depends on quantum interfer- 
ence and superposition and is not just a simulation of the QC experience. The 
flexibility is generated by a modular approach to all of the initializations, op- 
erators, gates, and measurements, which then can be readily used to describe 
the basic QC Teleportation [4], Grover's search [5,6] and Shor's factoring [7] 
algorithms. We also adopt a density matrix approach as an organizational 
framework for introducing fundamental Quantum Computing concepts in a 
manner that allows for more general treatments, such as handling the dynam- 
ics stipulated by realistic Hamiltonians and including environmental effects. 
That approach allows us to invoke many of the dynamical theories based on 
the time evolution of the density matrix. Since much of the code uses the 
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density matrix, we call it "QDENSITY," which stands for Quantum computing 
with a density matrix framework. However, the code also provides the tools to 
work directly with multi-qubit states as an alternative to the density matrix 
description. 

In section 2, we introduce one qubit state vectors and associated spin opera- 
tors, including rotations, and introduce commands from QDENSITY. The basic 
characteristics of the density matrix are then discussed in a pedagogic manner 
in section 3. Then in section 4, methods for handling multi-qubit operators, 
density matrices, and state vectors with commands from QDENSITY are pre- 
sented. Also in that section, we show how to take traces, subtraces and how 
to combine subtraces with general projection operators to simulate projective 
measurements in a multi-qubit context. The basic one, two and three qubit 
gates (Hadamard, CNOT, CPHASE, Toffoli, etc.) needed for the QC circuits 
are shown in section 5. The production of entangled states, such as the two- 
qubit Bell [8] states, the three-qubit GHZ [9] states, among others, [10] arc 
illustrated in both density matrix and state vector renditions in section 6. In 
sections 7-9, Teleportation, Grover's search, and Shor's factoring algorithms 
are outlined, with the detailed instructions relegated to associated notebooks. 
Sample application to the cluster or "one-way computing" model of QC is 
presented in section 10. Possible future applications of QDENSITY are given in 
the conclusion section 11. 

The basic logical gates used in the circuit model of Quantum Computing are 
presented in a way that allows ease of use and hence permits one to construct 
the unitary operators corresponding to well-know quantum algorithms. These 
algorithms are developed explicitly in the Mathematica notebooks as a demon- 
stration of the application of QDENSITY. A tutorial notebook (Tutorial. nb) 
available on our web site guides the user through the requisite manipulations. 
Many examples from QDENSITY, as defined in the package file Qdensity.m, are 
discussed throughout the text, which hopefully, with the tutorial notebook, 
will help the user to employ this tool. 

All these examples are instructive in two ways. One way is to learn how 
to handle QDENSITY for other future applications and generalizations, such as 
studying entanglement measures, examining the time evolution generated by 
experiment-based realistic Hamiltonians, error correction methods, and the 
role of the environment and its affect on coherence. Thus the main motivation 
for emphasizing a density matrix formulation is that the time evolution can be 
described, including effects of an environment, starting from realistic Hamil- 
tonians [11,12]. Therefore, QDENSITY provides an introduction to methods that 
can be generalized to an increasingly realistic description of a real quantum 
computer. 

Another instructive feature is to gain insight into how quantum superposition 
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and interference are used in QC to pose and to answer questions that would be 
inaccessible using a classical computer. Thus we can form an initial description 
of a quantum multi-qubit state and have it evolve by the action of carefully 
designed unitary operators. In that development, the prime characteristics of 
superposition and interference of probability amplitudes is cleverly applied in 
Quantum Computing to enhance the probability of getting the right answer 
to problems that would otherwise take an immense time to solve. 

In addition to applying QDENSITYto the usual quantum circuit model for QC, 
we have adapted it to examine the construction of cluster states and the steps 
needed to reproduce general one qubit and two qubit operations. These cluster 
model examples are included to open the door for future studies of the quite 
promising cluster model [13] or "one-way computing" approach for QC. 

We sought to simulate as large a system of qubits as possible, using new 
features of Mathematica. Of course, this code is a simulation of a quantum 
computer based on Mathematica code run on a classical computer. So it is 
natural that the simulation saturates memory for large qubit spaces; after all, 
if the QC algorithms always worked efficiently on a classical computer there 
would be no need for a quantum computer. 

Throughout the text, sample QDENSITY commands are presented in 
sections called "Usage." The reader should consult Tutorial. nb for 
more detailed guidance. 



2 ONE QUBIT SYSTEMS 



The state of a quantum system is described by a wave function which in gen- 
eral depends on the space or momentum coordinates of the particles and on 
time. In Dirac's representation independent notation, the state of a system is 
a vector in an abstract Hilbert space | ^(t) >, which depends on time, but in 
that form one makes no choice between the coordinate or momentum space 
representation. The transformation between the space and momentum repre- 
sentation is contained in a transformation bracket. The two representations 
are related by Fourier transformation, which is the way Quantum Mechanics 
builds localized wave packets. In this way, uncertainty principle limitations 
on our ability to measure coordinates and momenta simultaneously with arbi- 
trary precision are embedded into Quantum Mechanics (QM). This fact leads 
to operators, commutators, expectation values and, in the special cases when 
a physical attribute can be precisely determined, eigenvalue equations with 
Hermitian operators. That is the content of many quantum texts. Our pur- 
pose is now to see how to define a density matrix, to describe systems with 
two degrees of freedom as needed for quantum computing. 
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Spin, which is the most basic two- valued quantum attribute, is missing from a 
spatial description. This subtle degree of freedom, whose existence is deduced 
by analysis of the Stern-Gerlach experiment, is an additional Hilbert space 
vector feature. For example, for a single spin 1/2 system the wave function 
including both space and spin aspects is: 

V(r 1 ,t)\sm a >, (1) 

where | s m s > denotes a state that is simultaneously an eigenstate of the 
particle's total spin operator s 2 = s 2 + s 2 + s 2 ., and of its spin component 
operator s z . That is 

s 2 | sm s >= h 2 s(s + 1) | sm s > s z \ sm s >= hm s \ sm s > . (2) 

For a spin 1/2 system, we denote the spin up state as | sm s >— >| |, \ >=| >, 
and the spin down state as | sm s >— >| |, —\ >=| 1 >. 

We now arrive at the definition of a one qubit state as a superposition of the 
two states associated with the above and 1 bits: 

| ^ >= a | > +b | 1 >, (3) 

where a =< | \I/ > and b =< 1 | ^ > are complex probability amplitudes for 
finding the particle with spin up or spin down, respectively. The normalization 
of the state < ^ | ^ >= 1, yields | a \ 2 + | b \ 2 = 1. Note that the spatial 
aspects of the wave function are being suppressed; which corresponds to the 
particles being in a fixed location, such as at quantum dots. 2 

An essential point is that a QM system can exist in a superposition of these 
two bits; hence, the state is called a quantum-bit or "qubit." Although our 
discussion uses the notation of a system with spin, it should be noted that the 
same discussion applies to any two distinct states that can be associated with 
| > and | 1 >. Indeed, the following section on the Pauli spin operators is 
really a description of any system that has two recognizable states. 



2.0.1 Usage 

QDENSITY includes commands for qubit states as ket and bra vectors. For ex- 
ample, commands Ket[0], Bra[0],Ket[l], Bra[l], yield 



2 When these separated systems interact, one might need to restore the spatial 
aspects of the full wave function. 
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In[l] :=Ket[0] 
1 



Out[l] 







In[2] :=Ket[l] 
Out [2] 



V 



In[3] :=Bra[0] 
Out[3]:=(l 0) 



In[4] :=Bra[l] 
Out[4]:=(0 1) 



These are the computational basis states, i.e. eigenstates of the spin operator 
in the z-direction. 

States that are eigenstates of the spin operator in the x-direction are invoked 
by the commands 



In[5]:=BraX[0] 
which is equivalent to: 



In[6] := (Bra[0] + Bra[l])/V2 



Eigenstates of the spin operator in the y-direction are invoked similarly the 
commands BraY[0], BraY[l], etc. 
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2.1 The Pauli Spin Operator 



We use the case of a spin 1/2 particle to describe a quantum system with 
two discrete levels; the same description can be applied to any QM system 
with two distinct levels. The spin s operator is related to the three Pauli spin 
operators <J x ,<Jy, a z by 

(4) 

from which we see that 3 is an operator that describes the spin 1/2 system 
in units of |. Since spin is an observable, it is represented by a Hermitian 
operator, a' = 3. We also know that measurement of spin is subject to the 
uncertainty principle, which is reflected in the non-commuting operator prop- 
erties of spin and hence of the Pauli operators. For example, from the standard 
commutator property for any spin [s x , s y ] = ihs z , one deduces that the Pauli 
operators do not commute 

[a x , a y \ = 2ia z . (5) 
This holds for all cyclic components so we have a general form 3 

fa, Gj\ = 2ie ijk o- k . (6) 



An important property of the spin algebra is that the total spin commutes with 
any component of the spin operator [s 2 , s,] = for all %. The physical conse- 
quence is that one can simultaneously measure the total spin of the system 
and one component (usually s z ) of the spin. Only one component is a can- 
didate for simultaneous measurement because the commutator [s x , s y ] = ihs z 
is already an uncertainty principle constraint on the other components. As a 
result of the ability to measure s 2 and s z simultaneously, the allowed states of 
the spin 1/2 system are restricted to being spin- up and spin-down with respect 
to a specified fixed direction z, called the axis of quantization. States defined 
relative to that fixed axis are called "computational basis" states, in the QC 
jargon. The designation arises because as noted already one can identify spin- 
up with a state | >, which designates the digit or bit 0, and a spin-down 
state as | 1 >, which designates the digit or bit 1. 

The fact that there are just two states (up and down) also implies properties 
of the Pauli operators. We construct 4 the raising and lowering operators 

3 Here the Levi-Civita symbol is nonzero only for cyclic order of components 
ijk = xyz,yzx, zxy, for which = 1. For anti-cyclic order of components 
ijk = xzy, zyx,yxz = —1. It is otherwise zero. 

4 With the definition s± = s x ±is y , and using the original spin commutation rules, 
it follows that [s±,s 2 ] = =Fs±, which reveals that s± and hence also a± are raising 
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c± = &x ± ^Cy, and note that the raising and lower process is bounded 

ct+ | >= a- | 1 >= 0. (7) 

Hence, raising a two-valued state up twice or lowering it twice yields a null 
(unphysical) Hilbert space; this property tells us additional aspects of the 
Pauli operator. Since 

o±o± = [o x ± io y ) 2 = a 2 x - a 2 ± (o-^o-,, + (7^) = 0, (8) 

we deduce that = a 2 , and that the anti-commutator 

Wx,o- y } = a x a y + a y a x = 0. (9) 

The anti-commutation property is thus a direct consequence of the restriction 
to two levels. 

The spin 1/2 property is often expressed as: s 2 \ sm s >= h 2 s(s + 1) | sm s >= 
\% 2 | sm s >= | 2 (T 2 | sm s > . We have o 2 = 3 = a 2 x + a 2 + a 2 z = 2a 2 x + 1, where 
we use the above equality a 2 = a 2 , and from the z eigenvalue equation the 
property a 2 = 1, to deduce that 

o* x = oi = o* g = l. (10) 



Another useful property of the Pauli matrices is obtained by combining the 
above properties and commutator and anti-commutator into one expression 
for a given spin 1/2 system 

CTjfTj = 5ij + ie ijk <7 k , (11) 

where indices i, j, k take on the values x, y, z, and repeated indices are assumed 
to be summed over. For two general vectors, this becomes 

(a-A)(a-B) = A-B + i(AxB)-a. (12) 
For A — B — fj , a unit vector (a ■ fj) 2 = 1, which will be useful later. 

These operator properties can also be represented by the Pauli-spin matrices, 
where we identify the matrix elements by 

< s m s | o z | s m s > — > . (13) 



and lowering operators. The general result, including the limit on the total spin is 
s± | s m s >= y/s(s + 1) — m s (m s ± 1) | s m s ± 1 > . 
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Similarly for the x— and y— component spin operators 



< sm s | a x | s m s >- 




< sm s | <jy | sm s >- 



These are all Hermitian matrices <7j = aj 




Also, the matrix realization for the raising and lowering operators are: 



(14) 



(° 










a + = 


i 






:) 








i 2 





(15) 



Here 0+ = <7_. 

Note that these 2x2 Pauli matrices are traceless TV [a] = 0, unimodular of = 1 
and have unit determinant | det <jj |= 1. Along with the unit operator 



c = l 



1 



(16) 



the four Pauli operators form a basis for the expansion of any spin operator in 
the single qubit space. For example, we can express the rotation of a spin as 
such an expansion and later we shall introduce a density matrix for a single 
qubit in the form p = a + b-a = a + bn-a to describe an ensemble of particle 
spin directions as occurs in a beam of spin-1/2 particles. 



In QDENSITY, we denote the four matrices by <7j where i = is the unit matrix 
and % = 1,2,3 corresponds to the components x, y, and z. To produce the 
Pauli spin operators in QDENSITY, one can use either the Greek form or the 
expression s[i\. 



2.1.1 Usage 

QDENSITY includes commands for the Pauli operators. For example, there are 
three equivalent ways to invoke the Pauli a y matrix in QDENSITY: 



M7]:=<7 y 



Out [7] := 



-i 

1 
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In[8] :=s[2 
/ 

Out [8] : = 

V 



-i 

1 



The third way is to use the commands SigmaO, Sigmal,Sigma2, or Sigma3. 

Note that 



In [9] :=a 2 . KetY[0] - KetY[0] 



Out [9] := 



and 



In[10] := a 2 . KetY[l] + KetY[l] 



Out [10]: = 



confirm that KetY[0] and KetY[l] are indeed eigenstates of a y . Note that 
the • is used to take the dot product. 



2.2 Pauli Operators in Hilbert Space 



It is often convenient to express the above matrix properties in the form of 
operators in Hilbert space. For a general operator Q, using closure, we have 

n = ^2^2 \ n >< n \ n \ n' >< n' \ . (17) 

n n' 

For the four Pauli operators this yields: 



(7 = 

o x = 
o 2 = 
cr z = 



>< | + | 1 >< 1 | 
>< 1 | + | 1 >< | 

-i | >< 1 | +i | 1 >< 
>< I - I 1 >< 1 I . 



(18) 
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Taking matrix elements of these operators reproduces the above Pauli matri- 
ces. Also note we have the general trace property 

Tr | a >< b\=^2 <n \ a >< b | n >= < b | n X n \ a >=< b\ a >, 

n n 

(19) 

where | n > is a complete orthonormal (CON) basis 5 . Applying this trace rule 
to the above operator expressions confirms that Tr[a x ] = Tr[<r y ] = Tr[<r 2 ] = 0, 
and Tr[(j ] = 2. 

Another useful trace rule is Tr[ Q \ a >< b \ ] =< b \ Q \ a >. 



2. 3 Rotation of Spin 



Another way to view the above superposition, or qubit, state is that a state 
originally in the z direction | >, has been rotated to a new direction specified 
by a unit vector h = (n x ,n y ,n z ) = (sin 9 cos 0, sin 9 sin <f>, cos 9), as shown in 
Fig. 1. 





Fig. 1. Active rotation to a direction h (a); and active rotation around a vector fj 
(b). 

The rotated spin state 



n 



>= cos(9/2)e-^ 2 | > +sin(#/2)e+^ 2 | 1 >= 



cos(#/2)e-^/ 2 
[sm{9/2)e +t ^ 2 



is a normalized eigenstate of the operator 



(20) 



a ■ n = 



n x - m y 



n z 

n x + in y —n 



cos 9 sin9e %,j 
sin 9e l<t> — cos 9 



(21) 



We see that 



a ■ n n >=\ n > . 



(22) 



For a CON basis, we have closure Yl n \ n > <n\=l. 
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The half angles above reflect the so-called spinor nature of the QM state of 
a spin 1/2 particle in that a double full rotation is needed to bring the wave 
function back to its original form. 

These rotated spin states allow us to pick a different axis of quantization n for 
each particle in an ensemble of spin 1/2 systems. These states are normalized 
< n | n >= 1, but are not orthogonal < n' \ n >^ 1, when the n angles 9, (j) 
do not equal the n' angles 9', <fi'. 

Special cases of the above states with spin pointing in the directions ±x and 
±y are within phases: 




1 > ± | 1 > 
> ±i | 1 > 

71 



(23) 



Hilbert space versions are also shown above. 

Rotation can also be expressed as a rotation of an initial spin-up system about 
a rotation axis fj by an angle 7. Thus an operator R 1 = e - *^ "'^, acting as 



* >= i? 7 I > 



(24) 



can also rotate the spin system state to new direction. This rotation operator 
can be expanded in the form 



7 



7 



e — cos _ ao _ tsm _ a ■ i], 



(25) 



which follows from the property that (a-fj) 2 = fj-fj + i(fjxfj)-a = 1 . A special 
case of the state generated by this rotation i? 7 |0>isa7 = 7r/2 rotation 
about the fj — > y axis. Then the rotation operator is 



R-k/2 — e 4 " = cos — (7 — i sm — a y . 



(26) 



Introducing the Pauli matrices, this becomes 



Rtt/2 



1 -1 



V2 



(27) 
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This rotation about an axis again yields the same result; namely, 



Rn/2 \ 0> =^ 



i -i Wi = Wi 
l l V2 o y/2 i 



(2S 



Similar steps apply for a rotation about the x axis by 7 = 7r/2, which yields 
the earlier | ±y > states. 

From normalization of the rotated state < * | * >= < | R\R^ | > = < | > 
1, we see that the rotation is a unitary Rl f R 1 = 1 operator. 



2.3.1 Usage 

The trace of the Pauli operators is invoked in QDENSlTYby: 



InflJ^Tr^] 
Outjl] :=0 



Rotation about the X axis is represented by 



In[2]:=RotX[0] 
Out [2] : 



( Cos 


~e~ 

_2_ 




— iSin 


' e' 

_2_ 




I — ?Sin 


_2_ 


Cos 


~e~ 

2_ 





Commands for other directions are described in Tutorial. nb; see RotX[#], 
RotY[0], Rotqbit[vec,0]. 



2.4 One Qubit Projection 



For a one qubit system, it is simple to define operators that project on to the 
spin-up or spin-down states. These projection operators are: 



V =10x0 



V x =1 1 >< 1 



(29) 
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These are Hermitian operators, and by virtue of closure, sum to one J2 a =o,i 'Pa = 1- 
They can also be expressed in terms of the a z operator as 



l + a z 



v 1 



1 — (T, 



(30) 



or in matrix form 



Vn 



1 






1 



(31) 



One can also project to other directions. For example, projection of a qubit 
on to the ±x or ±y directions involves the projection operators 



V± x = | ±x > < ±x 
V± y = \ ±y >< ±y 



2 ' 
1 ± cr„ 



(32) 



2.4.1 Usage 

The above projections operators are invoked in QDENSlTYby: 



In[l]:=Po 



Out[l] := 



1 




In[2]:=P! 
/ 

Out [2] : = 

V 




1 



Projection operators using the x-basis are invoked by: 



In[3] := 
Out [3] := 



PXfO 
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A general operator ProjG[a,vec] to project into a direction stipulated by a 
unit three vector vec, with a=0 or 1, is also provided. See Tutorial. nb for more 
examples. These operators are useful for projective measurements. 



3 THE (SPIN) DENSITY MATRIX 



The above spin rotated wave functions can be used to obtain the expectation 
value of some relevant Hermitian operator Q = fl\ which represents a physical 
observable. Let us assume that a system is in a state labelled by a with a state 
vector | a >. In general the role of the label a could be to denote a spatial 
coordinate (or a momentum), if we were considering an ensemble of localized 
particles. For the spin density matrix, we use a to label the various spin 
directions n. 

The average or expectation value of a general observable Q is then < a \ Q \ a > 
This expectation value can be interpreted simply by invoking eigenstates of 
the operator Q 

Q | v >= cu u | v >, (33) 

where u v are real eigenvalues and | v > are the eigenstates, which are usually 
a complete orthonormal basis (CON). The physical meaning of the eigenvalue 
equation is that if the system is in the eigenstate | v >, there is no uncertainty 
Afl in determining the eigenvalue, e.g. 

(Afi) 2 =< v | Q 2 | v > - < v | Q | v > 2 = ul - col = 0. (34) 



Using the eigenstates | v >, we can now see the basic meaning of the expecta- 
tion value, which is a fundamental part of QM. The eigenstates form a CON 
basis. That means any function can be expanded in this basis and that the 
coefficients can be obtained by an overlap integral. For example, in general 
terms the completeness (C) allows the expansion 

l^>=E c -l z/> - ( 35 ) 

V 

The OrthoNormal (ON) aspect is < v \ v' >= 8 VV >. Thus 

< v' | ^ >= J2 c u < v' | v > = cv , (36) 
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reinserting this yields 

| ^ >=Y1 <u I ^ > \ u >= E I v >< v I ^ >= 1 1 ^ > • 

V V 

Thus we see completeness with orthonormality of the basis can be expressed 
in the closure form 

E \ v><v\=\ (37) 

V 

with I the unit operator in the Hilbert space. 

With closure (e.g. a CON basis), we can now see that the expectation value 
breaks in to a sum of the form 

< a \ Q \ a > = ^2^2 < a \ v >< v \ Q \ v' >< v' \ a > 
= ^2uj u < v | a >< a | v >= u u P^. 



Here =< u \ a >< a \ u >=\< u \ a >\ 2 \s the positive real probability of 
the state | a > being in the eigenstate | v >. Hence we see that the quantum 
average or expectation value is a sum of that probability times the associated 
eigenvalue u u over all possible values v. That is the characteristic of a quantum 
average. 

As the next step towards the spin density matrix, consider the case that we 
have an ensemble of such quantum systems. Each system is considered not 
to have quantum interference with the other members of the ensemble. That 
situation can be realized by the ensemble being located at separate sites with 
non-overlapping localized wave packets and also in the case of a low density 
beam, i.e. separated particles in the beam. This allows us to take a classical 
average over the ensemble. 

Suppose that the first member of the ensemble is produced in the state | a >, 
the next in | a' >, etc. The ensemble average is then a simple classical average 

< n >= J a>p « , (38) 

2-ta a 

where P a is the probability that a particular state a appears in the ensem- 
ble. Summing over all possible states of course yields X^ a P a = 1- The above 
expression is a combination of a classical ensemble average with the quantum 
mechanical expectation value. It contains the idea that each member of the 
ensemble interferes only with itself quantum mechanically and that the en- 
semble involves a simple classical average over the probability distribution of 
the ensemble. 
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We are close to introducing the density matrix. This is implemented by using 
closure and rearranging. Consider 



< a | Q | a > P a = < a \ m >< m \ Q \ m >< m \ a > P Q , 

a a mm i 

(39) 

where | m > denotes any CON basis. Now rearrange the above to 

< a | Q | a > P a = X! < m ' I a >< a \ m > Pa < m\Vt\m! > , 

a mm' a 

(40) 

and then define the density operator by 

p = | « >< « | P a (41) 

a 

and the associated density matrix in the CON basis | m > as < m | p | m' >= 
J2 a <m\a><a\m'> P a . We often refer to either the density operator 
or the density matrix simply as the "density matrix," albeit one acts in the 
Hilbert space and the other is an explicit matrix. The ensemble average can 
now be expressed as a ratio of traces 6 

<«>- W' (42) 

which entails the properties that 

Tr[p] = X! < m I P I m >= X! P« y j <a\m >< m \ a > 

m am 

= £ P a <a\a>= P« = l, (43) 

a a 

and 

Tr[pf2] = < m \ p\m >< m \ Q \ m > 

mm 1 

= X! X! P Q < a | m >< m \Vt\m >< m \ a > 

a mm' 

= ]TP a <a|ft|a>, (44) 

a 

which returns the original ensemble average expression. 



6 The trace TV is defined as the sum of the diagonal matrix elements of an operator, 
where a CON basis is used. 
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3. 1 Properties of the Density Matrix 



We have defined the density operator by a sum involving state labels a for 
the special case of a spin 1/2 system. The definition 

p = J2 \a><a\P a (45) 

a 

is however a general one, if we interpret a as the label for the possible charac- 
teristics of a state. Several important general properties of a density operator 
can now be delineated. The density matrix is Hermitian, hence its eigenvalues 
are real. The density matrix is also positive definite, which means that all of 
its eigenvalues are greater or equal to zero. This, together with the fact that 
the density matrix has unit trace, ensures that the eigenvalues are in the range 
[0,1]. 

To prove that the density matrix is positive definite, consider a basis | v > 
which diagonalizes the density operator so that 



< v | p | v > — n v (46) 
= ^P a < v | a >< a | v >= l< v I a >| 2 > 0- 

a a 

Here p v is the z/th eigenvalue of p and both parts of the final sum above are 
positive quantities. Hence all of the eigenvalues of the density matrix are > 
and the density matrix is thus positive definite. If one of the eigenvalues is 
one, all the others are zero. 

Another general property of the density matrix involves the special case of a 
pure state. If every member of the ensemble has the same quantum state, then 
only one a (call it ceo) appears and the density operator becomes p = \ «o >< «o 
The state | cto > is normalized to one and hence for a pure state p 2 = p. Us- 
ing a basis that diagonalizes p, this result tells us that the eigenvalues satisfy 
Pu(Pu — 1) = and hence for a pure state one density matrix eigenvalues is 1, 
with all others zero. 

In general, an ensemble does not have all of its members in the same state, 
but has a mixture of possibilities as reflected in the probability distribution 
P a . In general, as we show below, we have 

P 2 < A (47) 

with the equal sign holding for pure states. A simple way to understand this 
relationship is seen by transforming the density matrix to diagonal form, using 
its eigenstates to form a unitary matrix U p . We have U p pU^ p = po, where pn 
is diagonal using the eigenstates of p as the basis, e.g. < u \ p D \ u' >= p u 5 uu >. 
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Here p v again denotes the z/th eigenvalue of p. We already know that the sum 
of all these eigenvalue equals 1, that they are real and positive. Since every 
eigenvalue is limited to be less than or equal to 1, we have p? v < p u , for all v. 
Transforming that back to the original density matrix yields the result p 2 < p. 
Taking the trace of this result yields another test for the purity of the state 
Tr[p 2 ] < Tr[p] = I. Examples of how to use this measure of purity will be 
discussed later. 

3.1.1 Entropy and Fidelity 

As an indication of the rich variety of functionals of p that can be defined, let 
us examine the Von Neumann entropy and the fidelity. 

The Von Neumann entropy [14], S[p] = — Tr[p log 2 p], is a measure of the 
degree of disorder in the ensemble. Its basic properties are: S[p] = if p 
is a pure state, and S[p] = 1 for completely disordered states. See later for 
an application to the Bell, GHZ , & Werner states and also the Tutorialnb 
notebook for simple illustrative examples. 

It is often of interest to compare two different density matrices that are alter- 
nate descriptions of an ensemble of quantum systems. One simple measure of 
such differences is the fidelity. Consider for example, two pure states 

p =\ ijj >< ijj \ p =| tp >< $ |, (48) 

and the associated overlap obtained by a trace method 

Tr[ p p ] = Tr[ | V >< i> \ 4> >< $ | ] = | < i> \ i> > | 2 - ( 49 ) 

Clearly this overlap equals one if the states and associated density matrices 
are the same and thus serves as a measure to compare states. This procedure 
is generalized and applied to general density matrices. It is also written in a 
symmetric manner, with the general definition of fidelity being 

F[p,p]=Tr[ ] l^pp^~p] , (50) 

which has the property of reducing to Tr[p] = 1, for p — p. It also yields 
| < ip | tp > | in the pure state limit. 

3.1.2 Usage 

QDENSITY includes commands that produce the Purity and Entropy for a stip- 
ulated density matrix p, Purity[p], Entropy [p], and the Fidelity of one specified 
density matrix relative to another Fidelity [p±, p?\. 
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3.1.3 Composite Systems and Partial Trace 

For a composite system, such as colliding beams, or an ensemble of quantum 
systems each of which is prepared with a probability distribution, the defini- 
tion of a density matrix can be generalized to a product Hilbert space form 
involving systems of type A or B 

PAB = I a @ >< a @ I' ( 51 ) 

where P a ,/3 is the joint probability for finding the two systems with the at- 
tributes labelled by a and (3. For example, a could designate the possible 
directions n of one spin-1/2 system, while /3 labels the possible spin directions 
of another spin 1/2 system. One can always ask about the state of system 
A or B by summing over or tracing out the other system. For example the 
density matrix of system A is picked out of the general definition above by 
the following trace steps 

Pa = ^b[pab] 

= p «,/3 \ a><a\Tr B [\P><P\] 

= H(H P «,/3) I a >< a I 

a 

= J2 P a I ct >< a I . (52) 

a 

Here we use the product space | a/3 >^\ a >\ j3 > and we define the proba- 
bility for finding system A in situation a by 

P« = E P ^- (53) 

3 

This is a standard way to get an individual probability from a joint probability. 

It is easy to show that all of the other properties of a density matrix still hold 
true for a composite system case. It has unit trace, it is Hermitian with real 
eigenvalues, etc. 

See later for application of these general properties to multi-qubit systems. 
3.2 Comments about the Density Matrix 
3.2.1 Alternate Views of the Density Matrix 

In the prior discussion, the view was taken that the density matrix implements 
a classical average over an ensemble of many quantum systems, each member 



22 



of which interferes quantum mechanically only with itself. Another viewpoint, 
which is equally valid, is that a single quantum system is prepared, but the 
preparation of this single system is not pinned down. Instead all we know is 
that it is prepared in any one of the states labelled again by a generic state 
label a with a probability P Q . Despite the change in interpretation, or rather 
in application to a different situation, all of the properties and expressions pre- 
sented for the ensemble average hold true; only the meaning of the probability 
is altered. 

Another important point concerning the density matrix is that the ensem- 
ble average (or the average expected result for a single system prepared as 
described in the previous paragraph) can be used to obtain these averages 
for all observables Q. Hence in a sense the density matrix describes a sys- 
tem and the system's accessible observable quantities. It represents then an 
honest statement of what we can really know about a system. On the other 
hand, in Quantum Mechanics it is the wave function that tells all about a 
system. Clearly, since a density matrix is constructed as a weighted average 
over bilinear products of wave functions, the density matrix has less detailed 
information about a system that is contained in its wave function. Explicit 
examples of these general remarks will be given later. 

To some authors the fact that the density matrix has less content than the 
system's wave function, causes them to avoid use of the density matrix. Others 
find the density matrix description of accessible information as appealing. 



3.2.2 Classical Correlations and Entanglement 

The density matrix for composite systems can take many forms depending on 
how the systems are prepared. For example, if distinct systems A & B are 
independently produced and observed independently, then the density matrix 
is of product form p A B l— ► Pa®Pb, and the observables are also of product form 
^ab l— ► <8> &b- For such an uncorrected situation, the ensemble average 
factors 

/ n ^ ^[Pab&ab] Tr[p A Q A ] Tr[p B fi B ] 
K ^ >= Tr[ PAB ] = Tr M Tr[p B ] (54) 
as is expected for two separate uncorrelated experiments. This can also be 
expressed as having the joint probability factor P a ,/3 l— * P«P/3 the usual prob- 
ability rule for uncorrelated systems. 

Another possibility for the two systems is that they are prepared in a coordi- 
nated manner, with each possible situation assigned a probability based on the 
correlated preparation technique. For example, consider two colliding beams, 
A & B, made up of particles with the same spin. Assume the particles are pro- 
duced in matched pairs with common spin direction n. Also assume that the 
preparation of that pair in that shared direction is produced by design with a 
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classical probability distribution P T - t . Each pair has a density matrix pn <g> pn 
since they are produced separately, but their spin directions are correlated 
classically. The density matrix for this situation is then 

Pab = ^2Ph ph® ph- (55) 

h 

This is a "mixed state" which represents classically correlated preparation and 
hence any density matrix that takes on the above form can be reproduced by 
a setup using classically correlated preparations and does not represent the 
essence of Quantum Mechanics an entangled state. 

An entangle quantum state is described by a density matrix (or by its cor- 
responding state vectors) that is not and can not be transformed into the 
two classical forms above; namely, cast into a product or a mixed form. For 
example, a Bell state \(\ 01 > + | 10 >) has a density matrix 

p = -(| 01 >< 01 | + | 01 >< 10 | + | 10 >< 01 | + | 10 >< 10 | ) (56) 

that is not of simple product or mixed form. It is the prime example of an 
entangled state. 

The basic idea of decoherence can be described by considering the above case 
with time dependent coefficients 

p = -( ai {t) | 01 >< 01 | +a 2 (t) | 01 >< 10 | +a* 2 (t) | 10 >< 01 | +a 3 (t) | 10 >< 10 | ). 

2 (57) 
If the off-diagonal terms a 2 (t) vanish, by attenuation and/or via time averag- 
ing, then the above density matrix does reduce to the mixed or classical form, 
which is an illustration of how decoherence leads to a classical state. 



4 MULTI -QUBIT SYSTEMS 



The previous discussion which focused on describing a single qubit, can now 
be generalized to multiple qubits. Consider the product space of two qubits 
both in the up state and denote that product state as|00>=|0>|0>, 
which clearly generalizes to 

I <?1 <?2 >=| <?i >| <?2 >, (58) 

where q±, q 2 take on the values and 1. This product is called a tensor product 
and is symbolized as 

I <?i <?2 >=| <?i > ® I <?2 >, (59) 
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in[5]:- ketl23= ( (KetX[0] 8Ket [1] ) SKetX [1] ) 



ln[4]:= ketl2 = KetY[0] ®Ket[l] 



In[4]:= KetV[{0, 1, 1}] 



Out [4] 



Fig. 2. Simple examples of tensor products of two and three kets. 
which generalizes to n q qubits 



Q1Q2 ■■■ Qn q >= (\qi > ® I <?2 >) 



(60) 



In QDENSITY, the kets | >, | 1 > are invoked by the commands Ket[0] and 
Ket[l], as shown in Fig. 2, along with the kets | ±x >, and | ±y >. 

Also shown in that figure are the results for forming the tensor products of 
the kets for two and three qubits, as described next. 



4-1 Multi-Qubit Operators 



One can also build operators that act in the multi-qubit spin space described 
above. Instead of a single operator, we have a set of separate Pauli operators 
acting in each qubit space. They commute because they refer to separate, 
distinct quantum systems. Hence, we can form the tensor product of the n q 
spin operators which for two qubits has the following structure 



< ai ] <7j | bi >< a 2 | o~j | b 2 > = < a\a 2 \ <j\ a) \ b\b 2 > 

= < a x a 2 | af ] ® of \ b x b 2 > , (61) 

which defines what we mean by the tensor product of £g> of for two qubits. 
The generalization is immediate 

(of®of)®(of®of)--- . (62) 

The corresponding steps in QDENSITY are shown in Fig. 3, 

For large numbers of qubits a recursive method has been developed (see, Qden- 
sity.m and Tutorial.nb) , which involves specifying the "Length" L= n q of the 
qubit array and an array of length L that specifies the Pauli components 
k,l,- ■ ■. For example, if i — l,j — there is a o x in qubit 1 space and 
a unit operator o acting in qubit 2 space. The multi-qubit spin operator is 
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Fig. 3. Multi-qubit operators in QDENSITY 

called SP[L, {i, j, k, I, • • • }]. Examples in Fig. 3 include operator tensor prod- 
ucts generated directly using the ® notation. 

4-1.1 Usage 

QDENSITY includes a multiqubit spin operator SP[L,{a 1 ,a 2 v? a L}] built from 
L Pauli spin operators of components ai,a 2 ,..,QL. A sample construction is: 



In[l]:=SP[2,{2,3}] 



Out[l] := 



t -i \ 



i 

1 
-i 



V 



which is equivalent to the tensor product 
In[2] 



Out [2] : = 



a 2 a s 

t -i 

i 

1 
y -% 



The advantage of this command is that it can readily construct large space 
tensor products. 
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4-2 General Multi -Qubit Operators 



The production of n q spin space operators provides a complete basis for ex- 
pressing any operator. This remark is similar to, and indeed equivalent to, the 
statement that the n q ket product space is a CON basis. With that remark, 
we can expand any n q operator as 

« = E ® ° { 2 ® <7g> • • • <> = E Ca SP[n 9 , a] (63) 

a a 

where the sum is over all possible values of the array a : {a±, 03, • • • , a nq }. 
Here, the multi-qubit spin operator is denoted by SP[n g , a], which is the no- 
tation used in QDENSITY. The coefficient C a can be evaluated for any given fl 
from the overall trace 

C a = ^-Tr[fi.SP[n 9 ,a]]. (64) 

Because of the efficacy of Mathematica 5.2, the total trace can be evaluated 
rapidly. This set of coefficients characterizes the operator Q. 



4-2.1 Partial Traces 

The advantage of expanding a general operator in the Pauli operator basis 
is that partial traces can now be generated by manipulating the above coef- 
ficients. A partial trace involves tracing out parts of a system; for example, 
consider the partial trace over qubit two for a three qubit operator 

Tr 2 K (1) ® af ® af\ = 2S j0 af ] ® af . (65) 

Recall that Tr[<7j] is zero unless we have the unit matrix <To in which case the 
trace is two. Of course, one could trace out systems 2 and also 1, and then 

Tr 12 [aV ® af ® af] = 25 l0 28 j0 af . (66) 

The subscript on the Tr symbol indicates which qubit operators are being 
traced out. Note in this case the number of qubits in the result is reduced to 
n q — 2, where 2 is the length of the subscript array in TV12. Clearly, the trace 
reduces the rank 7 of the operator by the number of qubits traced out. 

Now we can apply these simple ideas to construct the partial trace of a general 
operator f2. Using the linearity of the trace 

Tr t [fl] = £ C a Tr t [a« ® a® ® a% . . . a M] , (67) 

a 

7 The rank is the number of qubits n q . 
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In[7] := p 123 = ( (PX[0] ®PX[1] ) ®PY[0] ) 








Pl = PTr[{2, 


3}, P123] 












P12 = PTr[{l}, pi 23 ] 
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Fig. 4. Taking partial traces with QDENSITY 



where the array t : {tit 2 • • • } indicates only those qubits that are to be traced 
out. For example, t : {25} indicates that only qubits 2 and 5 are traced out. 

The procedure for taking a partial trace of a general operator is to deter- 
mine the total coefficient C a for all of the array a : {aid 2 • • -a^}, except for 
the entries corresponding to the traced out qubit for which we need only the 
cij = part if say we trace out the jth qubit. From the resultant coefficients, 
we obtained a reduced set of coefficients, reduced by the number of trace outs. 
That reduced coefficient is then used to construct the reduced space opera- 
tor, with a multiplier of 2 included for each traced out qubit. This expansion, 
reduction, reconstruction procedure might seem complicated, but it has been 
implemented very efficiently using the power of Mathematica 5.2. See Qden- 
sity.m for the explicit construction procedure (which is rather compact). The 
command used in QDENSITY is PTr [ t, Q] where the trace out of the general 
operator fl is specified by the array t. Examples of the partial traces are in 
Fig. 4. 

4.2.2 Usage 

QDENSITY includes several commands for taking partial traces. One is PTr[{g 1 ,g 2 ,---?<2 , M}?^] ) 
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PTr[{l, 2, 3}, SP[5, {1, 3, 3, 0, 1}]] 

/ \ 


.0 , 

PTr[{3, 4, 5}, SP[5, {1, 3, 0, 0, 0}]] 

'0 8 \ 

-8 

8 

.0 -8 j 

Fig. 5. Partial Traces of multi-qubit Pauli operators 

where the array qi,q2,---,qM stipulates the space to be traced out. See Tuto- 
rial. nb and Fig. 4 for examples of these commands. 

4-3 Multi-Qubit Density Matrix 

The multi-qubit density matrix is our prime example of an operator that we 
examine in various ways, including taking partial traces. Just as in the prior 
discussion, a general density matrix can be expanded in a Pauli spin operator 
basis 

P = E SPKa], (68) 

a 

where the coefficient C Pa is real since the density matrix and the Pauli spin 
tensor product SP[n q , a] are Hermitian. Taking a partial trace follows the rules 
discussed earlier. Examples are presented in Figs. 5 and 6. 

In these examples, we give the case of three qubits reduced to two and then 
to one. The general expansion for these cases takes on a simple and physically 
meaningful form and therefore are worth examining. For one qubit, the above 
expansion is of the traditional form 

Pl = i[l + /W], (69) 

which involves the three numbers contained in the vector Pi, also know as 
the polarization of the ensemble. A 2 x 2 Hermitian density matrix has 4 
variables, which is reduced by one by the Tr[pi] = 1 normalization. Thus the 
polarization vector is a complete parametrization of a single qubit. For a pure 
state, the magnitude of the polarization vector is one; whereas, the general 
constraint p 2 < p implies that | Pi |< 1. A graph of that vector thus lies within 
a unit circle called the Bloch sphere. The physical meaning of Pi is that it is 
the average polarization of an ensemble, which is made clear by forming the 



Out [10]= 



In [11] := 



Out [11]= 
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ensemble average of the Pauli spin vector: 



Now consider two qubits. The Pauli basis is U{ ® crj, and hence the two qubit 
density matrix has the form 



p 12 = - [ 1 + Pi • a x <S> 1 + 1 <S> o 2 • P 2 + era <S> cr 2 jT itj \ 
= I[ i + p x . B 1 + P 2 ■ a 2 + a x ■ *T ■ a 2j }. 



This involves two polarization vectors, plus one 3x3 tensor polarization T 
8 which comes to 15 parameters as indeed is the correct number for a two 
qubit system 2 2 x 2 2 — 1 9 . The physical meaning is again an ensemble aver- 
age polarization vector for each qubit system, plus an ensemble average spin 
correlation tensor 



Tr[p 12 gi g 1 2 ] _ £ 
<(Xl> - Tr[ Pl2 ] =Pl ' 

<a2> ~ Tr[ Pl2 ] =P2 ' 



Tr[p i 2 o- lf (8) q- 2j ] _ 
Tr[p 12 ] 

(71) 



< 0"li(7 2j > = — r— 1 = 



To illustrate a partial trace, consider the trace over qubit 2 of the two qubit 
density matrix 

Tr 2 [ Pl2 ]= Pl = i[l + iWi], (72) 

where we see that a proper reduction to the single qubit space results. Exam- 
ples of the density matrix for the Bell states and their partial trace reduction 
to the single qubit operator are presented in Fig. 6. 



8 In the tensor term the sum extends only over the i,j = 1, 2, 3 components 

9 We see that the number of parameters in a n q qubit density matrix is thus 2 nq x 
2 n 9 — i = 2 2n 9 — 1 
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In[2] 


= Bellstate 


= (Ket[0] ®Ket[0] +Ket[l] ®Ket[l]) /Sqrt[2] 




pBell 


= Bellstate®Adj [Bellstate] 




Entropy [pBell] 


Out [2] 
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Out [4] 


= 








The Bell state has entropy 0, is a pure state. 


In [5] :- 


pBl = PTr[{2), pBell] 
pB2 = PTr[{l}, pBell] 
Entropy [pBl] 
Entropy [pB2] 




Out[5]= 
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Out [7] = 
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Out [8] = 
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...but, the entropy of each of the subsystems is 1 , they are in a completely disordered state! 



Fig. 6. Example from Tutorial.nb. 



4-4 Multi-Qubit States 



The procedure for building multi-qubit states follows a path similar to our 
discussion of operators. First we build the computational basis states, which 
are eigenstates of the operator <g> £g> • • •a^ lq \ These states are speci- 
fied by an array a : {a±, 0,2 ■ ■ ■ a 1lq } of length n q , where the entries are either 
one or zero. That collection of binary bits corresponds to a decimal number 
according to the usual rule a\a^ ■ ■ ■ a nq — > a\2 nq + g^™ 9-1 + a nq 2°. The corre- 
sponding product state | a\a^ • • • a n >= | a\ > <E> | a>2 > <8> • • • | On > can be 
constructed using the command KetV[{ai, a 2 , ..}]. Any single computational basis 
state consists of a column vector with all zeros except at the location counting 
down from the top corresponding to its decimal equivalent. Examples of the 
construction of multiqubit states in QDENSlTYare given in Fig. 2. This capa- 
bility allows one to use QDENSITY without invoking a density matrix approach, 
which is often desirable to reduce the space requirements imposed by a full 
density matrix description. 
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4-4-1 Usage 



QDENSITY includes commands for one qubit ket vectors in the computational 
and in the x- and y-basis Ket, KetX, KetY, and also multiqubit product 
states using the command KetVfvec]. Example of its use is 

In[l] :=KetV[{0,l,l}] 



1 




V / 



Out[l] 



which is equivalent to 



In[l] := (Ket[0] <g> Ket[l]) <g> Ket[l] 



Out[l] := 





1 






V / 



5 CIRCUITS & GATES 



Now that we can construct multi-qubit operators and take the partial trace, 
we are ready to examine the operators that correspond to logical gates for 
single and multi-qubit circuits. These gates form the basic operations that are 
part of the circuit model of QC. We will start with one qubit operators in a one 
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qubit circuit and then go on to one qubit operators acting on selected qubits 
within a multi-qubit circuit. Then two qubit operators in two and multi-qubit 
situations will be presented. 



5.1 One Qubit Gates 



5.1.1 NOT 

The basic operation NOT is simply represented by the a x matrix since a x | >=| 1 > 
and a x I 1 >=| > . 



5.1.2 The Hadamard 

For the Hadamard, we have the simple relation 



«-^75 Ml- (73 > 



which can also be understood as a rotation about the fj = axis by 7 = n 
since 

K — e 2 ' = cos — tTo — 1 sin — a • rj — > — 2 = — . (74) 

2 2 y 2 

The Hadamard plays an important role in QC by generating the qubit state 

from initial spin up or spin down states, i.e. 

tt|0>= |0> + |1:> W|l>= |0> - |1> . (75) 



Having a Hadamard act in a multi-qubit case, involves operators of the type 
7i <S> 1 ® 7Y, for which Hadamards act on qubits 1 and 3 only. The command 
for this kind of operator in QDENSITYis Had[n g ,Q] where n q is the total number 
of qubits and the array Q : (71,(72, ••• of length n q indicates which qubit is or 
is not acted on by a Hadamard. The rule used is if qi > 0, then the ith qubit 
is acted on by a Hadamard, whereas qj = designates that the jth qubit is 
acted on by a unit 2x2 operator. For example, Had[3, {1,0,1}] has a Hadamard 
acting on qubits 1 and 3 and a unit 2x2 acts on qubit 2, which is the case 
given above. To get a Hadamard acting on all qubits, include all qubits in 
Q, e.g., use Q={1,1,1,....}. Thus, an operator HALL[L]=Had[L, {1,1,1....}] is also 
implemented where the array of l's has length n q of all the qubits. Another 
QDENSITY command had[n g ,q] is for a Hadamard acting on one qubit q out 
of the full set of n q qubits. So QDENSITY facilitates the action of a one qubit 
operator in a multi-qubit environment. 
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NiceRotation = RotX [6i ] ®RotY [e 2 ] 

Cos [ ] Cos [£§-] - Cos [ ■%- ] Sin [ 2j- ] - i Cos [ Sj. ] sin [ ] i Sin [ ] sin [ 4f- ] 

Cos[-^-] Sin[-^-] Cos[-^-] Cos[-^-J - i Sin [ ^- ] Sin [ ^- ] - i Cos [ ^- ] Sin [ ^- ] 

-i Cos[4f-] Sin[-^-] i Sin [ ^- ] Sin [ Sg. ] cos [ ^- ] Cos [ ^- ] - Cos [ ^- ] Sin [ ^- ] 

-i Sin[-^-] Sin[-^-] - i Cos [ ^j- ] Sin [ ^j- ] Cos [ ^- ] Sin [ ^ ] Cos[^-] Cos[^-] 



Fig. 7. Example of multiqubit rotation using QDENSITY. 
5.1.3 Usage 



QDENSITY includes several Hadamard commands. For single qubit cases use 
either 7i or had [1,1]. For a Hadamard acting on a single qubit within a set 
of L qubits use had[L,q]; for a set of Hadamards acting on selected qubits 
use Had[L, { 0, 1, 0, 1 • ■ • } ], and for Hadamards acting on all L qubits use 
HALL[L]. These are demonstrated in the tutorial. 



5.1.4 Rotations 



One can use the rotation operator R to produce a state in any direction. 
A rotation about an axis fj is given in Eq. (25). For special cases, such as 
the x,y,z and 7 = n, the expanded form reduce to —ia x ,—icr y and — ia z , 
respectively. For a general choice of rotation, one can use the "MatrixExp" 
command directly, or use the spinor rotation matrix for rotation to angles 9, (p. 

For a multi-qubit circuit, the rotation operator for say qubits 1 and 3 can 
be constructed using the command R 11 (g> 1 <g> R l2 <E> 1 <8> • • • , with associated 
rotation axes. Examples from QDENSITY are given in Fig. 7. 



5.1.5 Usage 

Rotation commands for rorarions about the x-, y- or z- axis by an ankle 9 are 
included in QDENSITY: Rotx[0], RotY[0], Rotz[0] In addition, Rotqbit[v,t] builds 
the matrix corresponding to a rotation around a general axis axis v by an 
angle t. 



5.2 Two Qubit Gates 



To produce a quantum computer, which relies on quantum interference, one 
must create entangled states. Thus the basic step of two qubits interacting 
must be included. The interaction of two qubits can take many forms de- 
pending on the associated underlying dynamics. It is helpful in QC, to isolate 
certain classes of interactions that can be used as logical gates within a circuit 
model. 
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5.2.1 CNOT 



The most commonly used two-qubit gate is the controlled-not (CNOT ) gate. 
The logic of this gate is summarized by the expression 

CNOT | c,t >=\ c,t©c>, 

where c = 0, 1 is the control bit and t — 0, 1 is the target bit. In a circuit 
diagram the • indicates the control qubits and the © indicates the target 
qubit. 



The final state of the target is denoted as "t © c" where © addition is un- 
derstood to be modular base 2. Thus, the gate has the logical role of the 
following changes (control bit first) | 00 > i — > | 00 >; | 01 > i — > | 01 > 
; | 10 > i— > | 01 >;| 11 > i— > | 10 > . All of this can be simply stated using pro- 
jection and spin operators as 



CNOT[c,t] =| > c < | <g>I t + | 1 > c < 1 | ©a* 



(76) 



with c and t denoting the control and target qubit spaces. The CNOT , which 
is briefly expressed as CNOT = VqI + V\cr x , is called the controlled-not gate 
since NOT = a x . A matrix form for this operator acting in a two qubit space is 



CNOT = 



( I 0^ 
10 
1 
\0 1 

The rows & columns are ordered numerically as: 00, 01, 10, 11. 



(77) 
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The CNOT gate is used extensively in QC in a multi-qubit context. Therefore, 
QDENSITY gives a direct way to embed a two-qubit CNOT into a multi-qubit 
environment. The command is CNOT [n q ,c,t] where n q is the total number of 
qubits and c and t are the control and target qubits respectively as in the 
following examples: If the number of qubits is 6, and a CNOT acts with qubit 3 
as the control and 5 is the target, the operator (which is a 2 6 x 2 6 matrix) is 
invoked by the command CNOT[6, 3, 5]. The command CNOT[6, 5, 3] has 6 
qubits, with qubit 5 the control and 3 the target. The basic case in Eq.(77) is 
therefore just CNOT[2, 1, 2]. 



35 



5.2.2 CPHASE 



Other two qubit operators can now be readily generated. The CPHASE gate, 
which plays an important role in the cluster model of QC is simply a controlled 



CPHASE[c,t] =| > c < | <g>I t + | 1 > c < 1 



which in two qubit space has the matrix form 



CPHASE 



1 



1 







1 





o\ 




- 1 / 



(78) 



(79) 



The multi-qubit version is CPHASE[n ? , c, t], with the same rules as for the 
CNOT gate. 



5.2.3 Other Gates 

The generation of other gates, such as swap gates 10 , Controlled-io^ ( also 
known as a CROT gate) are now clearly extensions of the prior discussion. 

The swap gate swaps the content of two qubits. It can be decomposed in a 
chain of CNOT gates: 



l*i> 
1*2) 



1*2) 
l*i> 



Another example is CROT 

CROT[c,t] =| > c < I <g>I t + I 1 > c < 1 
which, in two qubit space, has the matrix form 



CROT = 



(l 








°1 





1 

















1 


1° 





-1 





(80) 



(81) 



(82) 



A multi-qubit version of CROT[n 9 , c, t], with the same rules as for the CNOT 
gate can easily be generated by a modification of Qdensity.m. 



10 



See the QDENSITY command Swap . 
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Indeed, the general case of a controlled-f2, where £1 is any one-qubit operator 
is now clearly 

CQ[c,t] =| > c < | <g>I t + | 1 > c < 1 | <8)fi*, (83) 
with corresponding extensions to matrix and multi-qubit renditions. 



5.2.4 Usage 

QDENSITY includes the following two-qubit operators, acting between qubits 
c(control) and t(target) imbedded in a system of L qubits: CNOT [L,c,t], CPHASE 
[L,c,t], ControlledX [L,c,t], ControlledY [L,c,t] and Swap [L,gi,g 2 ]- A generic two qubit 
operator within a multiqubit system involving operators Opl and Op2 is 
TwoOp[L,g 1 ,g 2 ,Opl,Op2]. 



5.3 Three Qubit Gates 

The above procedure can be generalized to three qubit operators. The most 
important three qubit gate is the Toffoli [l]gate, which has two control bits 
that determine if a unit or a NOT(<Xj,) operator acts on the third (target) bit. 
The projection operator version of the Toffoli is simply 

Toffoli = V <g> V <g> 1 + V <g> V x <g> 1 + V x <g> V <g> 1 + V x <g> V x <g> a x , (84) 

which states that the third qubit is flipped only if the first two(control) qubits 
are both 1. 

For a multi-qubit system the QDENSITY command Toffoli [n q ,q±, g 2 , 93] returns 
the Toffoli operator with qi and q 2 as control qubits, and q 3 as the target qubit 
within the full set of n q qubits. 

The Toffoli gate can be specialized or reduced to lower gates and is a universal 
gate. 



5.3.1 Usage 

QDENSITY includes a generic three qubit operator within a multiqubit system 
involving operators Opl, Op2, and Op3 : ThreeOp[L,gi, q 2 ,q3, Opl, Op2,Op3]. 
The Toffoli gate is a special case and is invoked by the command Toffoli [L,c,c,t], 
where c, c, t specifies the two control and the one target qubit out of the set 
of L qubits. 
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Let's produce a uniform superposition of 4 qubits: 
p^ll = HALL [4] . Pi . Adj [HALL [4] ] 



6 16 16 16 16 16 16 16 15 15 16 16 16 16 16 16 



6 16 16 16 16 16 15 15 15 15 16 16 16 16 16 16 



6 16 16 16 16 16 15 15 15 15 16 16 16 16 16 16 



6 16 16 16 16 16 15 15 15 15 16 16 16 16 16 16 



Fig. 8. Construction of a uniform four qubit state. 
6 SPECIAL STATES 

As a prelude to discussing QC algorithms, it is useful to examine how to pro- 
duce several states that are part of the initialization of a quantum computer. 
These are the uniform superposition, the two-qubit Bell, [8] three-qubit GHZ, 
[9] and Werner [10] states. 

6.1 Uniform superposition 

In many QC processes, the first step is to produce an initial state of the n q 
qubits that is a uniform superposition of all of its possible computational basis 
states. It is the initialization of this superposition that allows a QC to address 
lots of questions simultaneously and is often referred to as the "massively 
parallel" feature of a quantum computer. 

The steps start with a state of all spin- up | 0000 • • • >, then every qubit is acted 
on by a Hadamard Tt <g) Tt <g) Tt ® . . . , which is done by the QDENSITY command 
HALL[n ? ]. Thus each up state is replaced by Tt | >= ^ 0> J^ 1> , and we have the 
uniform superposition 

n 2 n 9-l 

| >= HALLK] | >= — - £ | x >, (85) 

x=0 

where x is the decimal equivalent to all possible binary numbers of length n g . 

An example of this process, including the associated density matrices, is in 
Fig. 8. 
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6.2 Bell States 



The singlet and triplet states are familiar in QM as the total spin zero and one 
states, with zero spin projection (M = 0). They are the basic states that enter 
into the EPR discussion and they are characterized by their "entanglement." 
Bell introduced two more combinations, so that for two spin 1/2 systems the 
four Bell states are: 



B 00 > = -^= | 00 > + | 11 > 
v2 

#oi > = -?= | 01 > + | 10 > 
v2 

Bio > = ^t= | 00 > - | 11 > 



V2 
1 

71 



Bii > = 4= | 01 > - | 10 >, (86) 



or m one line | B ab >= 75 \ 0b > +(-l) a | lb >, where q is the NOT[q] 
operation. 

A circuit that produces these states, starting from the state | ab > (a, b = 1, 0) 
consists of a Hadamard on qubit one, followed by a CNOT . The QDENSITY version 
is thus: 

B [a_, b_] :=CNOT[2,l,2]Had[2,{l,0}](Ket[a]®Ket[b]). 

The density matrix version involves defining the unitary transformation U = 
CNOT[2, 1, 2].Had[2, 1, 0] and an initial density matrix p l ah =| ab >< ab |, 
then evolving to the density matrix for each of the Bell states 

PT = U- P l b -Ul (87) 

In Fig. 9 part of this process taken from Tutorial. nb is shown. The tutorial 
includes a demonstration that the Bell states have zero polarization (as is 
obvious from their definition), and a simple diagonal form for the associated 
tensor polarization T . 

Another useful application shown in the tutorial is that taking the partial 
traces of the Bell state density matrices, yield non pure single qubit density 
matrices and that the associated von Neumann entropy defined by S[p] = — Tr[p log 
is zero for the Bell states, but 1 for the single qubit density matrices p\ and 
P2- Thus each qubit is in a more chaotic state, which physically means they 
have zero average polarization. 11 This property is an indication of the entan- 
glement of the Bell states. 
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Since many different state vectors can yield a net zero average polarization, it is 
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The Bell state operator of qubits i and j can be decomposed into the action of a Hadamard 
on qubit i and a CNOT gate on qubits i and j. 

In[51]:= Bellop[qi_, qj_] :=CNOT[2, qi, qj].had[2, qi] 

The density matrix for any of the four Bell states 

m[52]:- Pb o : = Bellop [ 1 , 2 ] . p 00 . Ad j [Bellop [ 1 , 2]]; 

Pboi := Bellop[l, 2] .p i .Adj [Bellop [1, 2]]; 

Pbio := Bellop[l, 2] .p 10 .Adj [Bellop [1, 2]]; 

Pen := Bellop[l, 2] .pu .Adj [Bellop [1, 2]]; 

Checks two Bell states 



In [56]:= PfiOO 

Pboi 

(i o o $\ 

Out[56]= 0000 

It \ i 
i o o o o\ 

i- \ 

Out[57]- , , 

° T T 
v , 

Fig. 9. Example from Teleportation.nb 

6.3 GHZ States 



Three qubit states that are similar in spirit to the Bell states, were introduced 
by Greenberger, M. A. Home, and A. Zeilinger [9]. The basic GHZ state is 



* >=| GHZ >= — ( | 000 > + | 111 >), 
v2 



(88) 



which may be written: 



GHZ >= U GHZ I 000 > 



(89) 



with U GHZ = CNOT[3,l,2].CNOT[3,l,3].had[3,l] which corresponds to the 
following circuit: 

|0> 



|0> 
|0> 





H 


1 


i 1 


i — 




r 








) 

€ 


3— 



A complete set of eight GHZ states can be produced by the step 



U, 



GHZ 



abc >=| GHZ a6c >= -^=(| 06c > +(-l) a | lbc >). (90) 

v2 



clear that the density matrix stores less information than in a state vextor, albeit 
realistic statistical, information. 
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The GHZ state may be prepared by acting on a state with three qLbits all iritialized to with 
the followi ng operator 



In[5S]:- SSZcp - CHQT[3 ( L, 3] .CHQT[3, 2].had[3, L] ; 

Then the density matrix forthis GHZstate is: 



In[59] :- pHSZ - SBZop. ( (T'o&T'o) SPoJ .Adj [HSZop] 



□ □ □ 



□ □□□□□□ 



□ □ □ 
□ □ □ 



Fig. 10. Construction of a GHZ state. Example from Tutorial.nb 

For all eight of these three qubit states, the associated density matrix can be 
formed H GHZ afec >< GHZ afjc | and are seen in Tutorial.nb to have a 
simple structure. Taking partial traces to generate the two qubit pn,Pi3,P23 
and single qubit p±,p2,P3 density matrices, we see that for these states every 
qubit has zero polarization and a simple structure for the pair and three 
qubit correlation functions. In addition, the entropy of these GHZ set of states 
is zero and the sub-entropies of the qubit pairs and single qubits are all 1, 
corresponding to maximum disorder. 



A sample GHZ realization in QDENSlTYis given in Fig. 10 



6.4 Werner States 



Another set of states were proposed by Werner [10]. They are defined in terms 
of a density matrix: 

Pw = Aps + (1 - A)p u <g> p u , (91) 

where < A < 1 is a real parameter. The limit A = 1 yields a completely 
entangled Bell state density matrix Pb =\ B a i, >< B a b |, whereas lower values 
of A reduce the entanglement. The A = limit gives a simple two qubit prod- 
uct p u ® p u , where p u = \ is the density matrix for a single qubit with zero 
polarization, i.e. it corresponds to a chaotic limit for each qubit. Therefore, 
the parameter A can alter the original fully entangled Bell state by intro- 
ducing chaos or noise. Therefore, the Werner state is called a state of noisy 
entanglement. 

The entropy of a two qubit Werner state as a function of A ranges from two 
12 for A = 0, to zero for A = 1. The entropy of the single qubit state is 1. See 
Tutorial.nb for a sample Werner QDENSITY realization. 

12 This corresponds to an entropy per qubit of 1. 
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7 TELEPORTATION 



To understand QC teleportation, let us first consider classical teleportation, 
which entails only classical laws. For example, suppose Alice measures a vase 
using a laser beam to measure its dimensions, shape, color and decoration. 
She therefore has a full binary description of the vase in 3-D. Bob has all the 
material to make another vase and, upon receiving the file that Alice sends 
him by computer, is able to make an exact copy of the vase. There are only 
local operations (LO) (measuring and sending by Alice and reconstruction by 
Bob) and a classical communication (CC); so this called a LOCC process. How 
does this differ from teleportation using Quantum Physics? 

In Quantum Mechanics, measurement affects the sample; as a result, after 
collecting the information to send to Bob, the original sample is no longer in 
its original state. In the classical case, one ends up with two identical copies. In 
the QC case, Bob has the only extant system. Another difference is that in the 
QC case, Alice and Bob share an entangled state, say an EPR or a Bell state, 
Alice entangles the original system with one member of the pair, then measures 
and by LOCC sends Bob her result. By virtue of the shared entanglement, 
information is shared by the LOCC and by the Quantum effect of sharing the 
entangled state. Some information is transmitted by a "Quantum channel." 
Therefore, the information sent by computer is less than needed in the classical 
case, because it is supplemented by the Quantum transfer of entanglement. 
The strange nature of Quantum transportation is thus no stranger than the 
EPR/Bell effect, which has been affirmed experimentally. 

To understand these general remarks, let us use QDENSlTYto examine three 
cases. 



7.1 One Qubit Teleportation 

Suppose Alice has one qubit qi in an unknown state | \1/ >= a | > +ai | 1 >, 
with an associated spin density matrix p =| ^ >< \I> | . In QDENSITY, such 
a state is generated randomly. Bob and Alice share a two qubit entangled 
state, which we take as one of the Bell states B q2m . This Bell state could be 
provided by an outside EPR purveyor, but for convenience let us assume that 
Bob produces the entangled pair and sends one member of the pair | q 2 > 
to Alice, as shown in Fig. 11. Alice then entangles her state | ^ >, using 
the inverse of the steps that Bob employed to produce the entangled pair, 
and then she measures the state of her | > <g> | q 2 >, which yields a single 
number between zero and three (or one of the binary pairs 00; 01; 10; 11). Alice 
transmits by CC that number to Bob, who then knows what to do to his qubit 
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q% in order to slip it over to the state | ^ >, without making any measurements 
on it to avoid its destruction. In the end, Alice no longer has a qubit in the 
original state, but by LOCC and shared entanglement, she is happy to know 
that Bob has such a qubit, albeit not made of the original material. The only 
material transmitted is the single member of the entangled pair. 



l«2> = 10) 

l?3> = 10} 



H 



H 



0{M) 



Fig. 11. Schematic circuit for one qubit quantum teleportation. 

In QDENSITY, the following detailed steps are shown in the file Teleportation. nb: 
construction of an initial random qubit for Alice, Bob's entanglement process, 
Alice's entanglement, measurement and CC and Bob's consequent actions. 
This process is described using the density matrix language. A rendition using 
the quantum states directly is easily generated. 



7.2 Two Qubit Teleportation 



7.2.1 Two EPR Teleportation 

A similar process can be invoked to teleport an initially unknown state of two 
qubits | ^ >= X)j j=o i a>ij \ QuQ2j >■ 111 the special case that the unknown 
state is one of the Bell states, it can be transported using the procedure 
shown in Fig. 12. Bob now prepares a three qubit entangled GHZ state using 
Hadamards and CNOTs and then sends one of the three qubits q$ over to Alice, 
who entangles that one with her original qi , q 2 qubits | \I/ >i 2 <8> | q$ > and then 
measures the state of her three qubits, with the result of a number between zero 
and seven ( e.g. one of the binary results 000; 001; 010, 011; 100; 101; 110; 111). 
She transmits that decimal number to Bob by CC ( a phone call say), who 
then knows what to do to put his two qubits q±, q% into the original state 
| \& > . Again all the steps are presented in detail in QDENSITY in the file 
Teleportation. nb. 



7.2.2 General Two Qubit Teleportation 

The two qubits can be in a more general state than in the above discussion 
which was restricted to being one of the Bell states. In this case, Bob needs 
to entangle four qubits by a chain of Hadamard and CNOT gates as shown 
in Fig. 13. He then sends two qubits q 3 , g 4 over to Alice who entangles her 
two with them. That is, she entangles | \& > 12 <8> | (fe > <8> | q± >, then 
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Fig. 12. EPR teleportation using a GHZ entangled state. 

measures them with a decimal number result that is between zero and 14 or a 
binary measurement of: 0000; 0001; 0010; 0011; • • • ; 1111. With that number, 
Bob knows what to do and places his two qubits into the original state | \& >. 
Again see Teleportation. nb for the detailed layout. 
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Fig. 13. Schematic circuit for two qubit quantum teleportation. 



8 GROVER'S SEARCH 



Assume you have a set of items and you want to find a particular one in the 
set which has been marked beforehand. Let us further restrict the problem 
by saying that you are only allowed to ask yes/no questions, e.g. "Is this 
item the item?" In a disordered database with N items and one marked item 
that problem would require on the order of N trials to find the marked item. 
Quantum mechanics allows the states of a system to exist in superpositions 
and thus in many cases permits one to parallelize the process in some sense. 
Grover [5] proposed an algorithm that lowers the number of trials needed to 
0(\/~N) by making clever use of interference and superposition. He based his 
ideas on an analogy to the multiple slit problem [6]. 
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Step 5: Bob state is computed 

Bob's two qubit state, made of qubits 5 and 6, is obtained by taking the partial trace of the full final state after 
Alice's measurement (here is written in general depending on the four values of alice measurement) 

ln[28]:= | Bob[a_, b_, o_, d_] : = Chop [ PTr [ ( 1 , 2, 3, 4} , 2 " 4 * P,, b , c , a . p E .P,, b , c , a ] ] ; 

Alice measurement is produced at random, the result of her measurement is P1 ,P2,P3,P4 



ln[29]:= 



PI = Random [Integer, {0, 1}]; 

P2 = Random [Integer, {0, 1}]; 

P3 = Random [Integer, {0, 1}]; 

P4 = Random [Integer, {0, 1}]; 
{PI, P2, P3, P4} 



Fig. 14. Example from Teleportation.nb 

8.1 The Oracle 

Grover's search algorithm relies on the use of an Oracle. The idea underlying 
the Oracle is the following: the Oracle is a function that can recognize a 
solution to a problem although it may not know how to solve the problem. 
In that sense a classical example could be a lock and a key, the problem of 
finding the proper key that would open a lock out of a bunch of keys illustrates 
the role of the lock as an Oracle: you select a key and you try it, the Oracle 
would tell you whether that was the correct key or not, but you cannot ask 
the lock to single out the key from the bunch. The essential difference between 
a classical and a quantum Oracle is that in the quantum case the Oracle can 
act on an input which is a superposition of all states. In our example that 
would mean that we can try a superposition of all the keys at the same time. 

The above description of the role of an Oracle takes on the explicit form 

ORACLE \x> N \y >i = | x > N \ y © f(x) > u (92) 



which involves an N qubit and a single qubit product space and a specified 
function f(x). The matrix form of the Oracle is 



< x \< y' | ORACLE | x >| y >= 5 X , X <y' \ y@ f(x) > . (93) 



Examples of this Oracle matrix for single and double marked items are given 
in detail in the Grover.nb notebook, where the "inversion about the mean" 
process is also presented and explained in detail. 
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8.2 One marked item 



A schematic description of the searching process is: 
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The basic steps, which are detailed in Grover.nb, consist of applying a Hadamard 
to all qubits in the register while setting the Oracle to recognize the marked 
item. Then we need to construct the Grover operator and apply it a num- 
ber of times. Finally we measure the output register which then will give the 
marked item with high probability. The probability of success depends on the 
number of times the Grover operator has been applied, as can be studied in 
the notebook. 

Grover.nb contains two examples, the first one has only one marked item in 
the full database. The size of the database, given by the number of qubits, can 
be varied, together with the number of times the Grover operator is applied. 



8.3 Two marked items 

The second example includes two marked items in the database (and may be 
generalized to the case of M marked items). Of course, one needs to enlarge 
the register and Oracle space. 



9 SHOR'S ALGORITHM 

Shor's factoring algorithm is the most celebrated quantum algorithm, partly 
due to its powerful use of quantum superposition and interference to tackle 
one of the problems upon which most of our secure bank transactions rely. The 
factoring algorithm, described in detail in the notebook Shor.nb, essentially 
first relates the problem of factoring a number, N, made up of two prime 
numbers, N=p x q, to the problem of finding the period of a particular function 
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mod (N), being x a coprime 13 to N smaller than N. Then finding the 
period is related to the computation of the Quantum Fourier Transform (QFT) 
(analog to the discrete Fourier transform), for which a very effective quantum 
algorithm exists. 

Schematically, the procedure is the following: build two registers, 1 and 2, with 
rii and n 2 qubits each, initially set to |0 >. Then Hadamard the first register, 
the state of the full system then reads: 

Y 2 n i-l 

$ = ^ E I <>i ® I > 2 . (94) 

z i=0 

Then we take a number x smaller than N and coprime with it and load the 
second register with the function x l mod (N), giving: 

Y 2™i-l 

* = — I i >i ® I xi mod N >2 • ( 95 ) 
2ni i=0 

At this point, a measurement is performed on the second register and then one 
applies the QFT to the first register. From that value measured in the first 
register, one is able, with a certain probability, to factor the original number 
N. 

A detailed study of the performance of the algorithm, e.g. analysis of proba- 
bilities of success depending on the size of register considered experimentally 
can be done within the notebook. A thorough theoretical description of the 
algorithm can be found in Refs. [7,15]. 

Shor.nb contains four slightly different approaches to the algorithm, mainly 
differing in the procedure used to compute the QFT. The most constructive 
case and also the one appropriate when studying possible noise effects on parts 
of the circuit, is the density matrix one. There the QFT (see also QFT.nb) is 
obtained using unitary operators in the same way as occurs experimentally. 
That full QFT treatment implies that the QM method is not practicable for 
heavy computing simulations, say when number of qubits goes above 10. 

Then there are three other cases, two of them using the state vector expression 
for the QFT with and without explicit construction of both registers. 

Finally the last example, makes use of the already implemented discrete Fourier 
transform in Mathematica. This example is thus useful when emphasis is on 

13 Two positive integers a and b are said to be coprime if they have no common 
factor other than 1, that is, their greatest common divisor is 1. 
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Step 6: Measure register 1, instead: study all possible results and their probabilites 

We store the possible measurements and their probabilities in variable MeasProb: 

ln[55]:= | MeaeProb = 

ln[56]:= I For [num = 0, num < 2 n l , £3 = Chop [Abs [QFTVec [ [num + 1] ] [ [1] ] ] *2] ; 
I MeasProb = Append [MeasProb, {num, Q}]; nura++] 

Plots of the possible measurements and their probabilities are quite illustrative: 




Plotl = ListPlot [MeasProb, Plot Joined -» True, 

AxesLabel -» {Measured number, Probability}, Image Size -» 400]; 
ListPlot [MeasProb, Plot Joined -» True , AxesLabel-* {Measurednumber, Probability}, 
PlotRange -* {0, 0.30}, IraageSize -» 400] ; 



Probability 
0.00025 



100 200 300 400 500 



Measured number 



Fig. 15. Example from Shor.nb 

studying larger numbers to check probabilities of success, without concern for 
the actual quantum mechanical way of building the QFT. In Fig. 15 a snapshot 
from Shor.nb is shown. 



10 CLUSTER MODEL 



An alternative to the circuit model for QC has been suggested in a series of 
recent papers [13]. The basic idea of this approach is to start with an initial 
state that is highly entangled by virtue of two-qubit CPHASE operations be- 
tween nearest neighbors in a cluster of qubits. The CPHASE operations could 
be generated somehow by Ising model spin-spin interactions. Once an appro- 
priate cluster is designed, then a carefully planned set of single qubit mea- 
surements are made in various directions. The results of those measurements 
are passed on by classical communications, until one reaches a final qubit, or 
set of qubits, from which a result can be deduced once a local correction in- 
volving Pauli operators and the binary results of the measurement is invoked. 
This method is being developed, with the procedures for general algorithms 
still being formulated. It is however novel and promising, especially since it 
involves single qubit measurements, can generate gates without use of mag- 
netic fields to rotate spins, and holds forth the promise of error stability. It 
does require however a large increase in the number of qubits. The fact that 
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l+> l+> 




Fig. 16. Cluster model of CNOTgate; see: Cluster.nb and [13]. 

measurement collapse out qubits and essentially destroys the initial state is 
why this approach is called "one-way computing." Of course, one can recon- 
struct the initial state and try again, which is typical of QC. One can also 
reuse qubits after measurement. 

It turns out that the modular nature of QDENSITY allows one to use it to demon- 
strate and test some of the basic ideas of the cluster model of QC. The first 
illustration is the basic transport process involving two qubits. The two qubits 
are initially in a | ^ > <S> | + > state, where | ^ > is a general unknown one 
qubit state. Then a CPHASE operator acts between the two qubits, followed by 
a measurement of qubit one in the | +x >=| + > direction with a result of 
either a = zero or one. The second qubit proves to be in the state a% | ^ >. It 
is simple to confirm this algebraically. It forms the basic building block of the 
cluster model. It is illustrated in QDENSITY using a density matrix example. 

Another example of using QDENSlTYfor cluster model studies is the simplest 
CNOTgate, whose cluster model implementation from Ref. [13] is shown in Fig. 
16. 

This example is worked out in detail in Cluster.nb. 



11 CONCLUSION 



This simulation affords opportunities for many applications and extensions. 
The basic operations and manipulations are formulated in a modular manner 
and, as the illustrations and tutorial demonstrate, one can formulate and an- 
swer many important QC questions. Application to dynamical theories based 
on master equations, including environmental effects, are one challenge. In- 
voking and testing measures for entanglement and probing the role of noisy 
entanglement, imperfect gates, and of error correction protocols are other po- 
tential applications. Extending the study to general types of measurements 
and to cluster model cases is also of considerable interest. Finally, the descrip- 
tion of real experimental situations by suitable Hamiltonians and studying the 
stability of QC algorithms could be an important role for future study using 
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The structure of the Cluster is essentially described in ClusterLinks, which 
contains the pairs of CPHASE interacting qubits: 



NCS = nq 

ClusterLinks = {{1, 2}, {2, 3}, {2, 4}}; 
Nlinks = Length[ClusterLinks] 

The 4 qubit initial state before building in the correlations is: 



InitialState = ( (ketl ®KetX [0] ) ®KetX[0] ) ®ket4 
pinitial = InitialState® Adj [InitialState] ; 
Chop[Tr [InitialState ® Ad j [InitialState] ] ] 



Fig. 17. Example from Cluster. nb showing preparation of the cluster state. 
QDENSITY. 
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Single Qubit Tools 


P. 


PX[B] 


?¥[■] 




Ket [■] 


Br a [ b] 


KefcX [ ■] 


DTcLA [ BJ 


Ke tY [ ■ ] 


BraY [ sj 


Ketl[«, ■] 


Bral[a, b] 


ProjG[ Bj {■}] 


Tr[B] 


RamdomQubitl 


Adj[B] 


RandoinKetl 


RandoinBral 


Rotqbit[{-.}, ■] 


RotX[B] 


RotY [ ■] 


RotZ[B] 


UoiiuiL [ ■.. ■] 


Fidelity[B, a] 


Purity [■] 


Entropy [b] 






MultiQubit Tools 


KetV[{«}] 


SP[B, {B}] 


PTr[{a}, ■] 


(■«>■) 


MeasureKet[a, {■}, {■}] 


RandomQubit2 


RaridomQuMtNlB] 


MultiQProj [ {a}, {b}] 


Common Gates 


ClfOT[B, l, ■] 


CPHASE[B, a, b] 


Had[B, {B}] 


CPHIlSEV[B, a, {B}] 


ControlledI[[a, u, a] 


ControlledTlBj m, a] 


had [a, ■] 


HftLL [ a] 


Toff old. [a, a, ■] 


Swap [a, a, b] 







Fig. A.l. The QDENSITY palette with the most useful commands 



A The QDENSITY palette 



The package comes together with a Mathematica Palette containing most of 
the commands in the package and also some symbols that may be useful. A 
snapshot of the palette is in Fig. A.l. 
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